Mixture- like behavior near a liquid-liquid phase transition in simulations of 

supercooled water 
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In simulations of a water- like model (ST2) that exhibits a liquid- liquid phase transition, we test for 
the occurrence of a thermodynamic region in which the liquid can be modelled as a two-component 
mixture. We assign each molecule to one of two species based on the distance to its fifth-nearest 
neighbor, and evaluate the concentration of each species over a wide range of temperature and 
density. Our concentration data compare well with mixture-model predictions in a region between 
the liquid-liquid critical temperature and the temperature of maximum density. Fits of the model 
to the data in this region yield accurate estimates for the location of the critical point. We also 
show that the liquid outside the region of density anomalies is poorly modelled as a simple mixture. 

PACS numbers: 64.30.-t, 64.70.Ja, 64.60.De 



The possibility that a liquid-liquid phase transition 
(LLPT) occurs in supercooled water and other tetrahe- 
dral liquids (e.g. silicon) continues to be a subject of 
investigation and debate [H [!]• In the LLPT proposed 
for water, two phases, a low-density liquid (LDL) and 
high-density liquid (HDL), become distinct below a crit- 
ical temperature Tc located in the supercooled regime. 
While a LLPT for water has yet to be confirmed exper- 
imentally, simulation studies have identified unambigu- 
ous LLPTs in several model tetrahedral liquids, includ- 
ing water |^], silicon [3], and tetrahedrally coordinated 
colloids @. 

Long before any discussion of LLPTs, there were recur- 
ring proposals that the thermodynamic anomalies of wa- 
ter, such as the density maximum, could be understood if 
the liquid is modelled as a mixture of two "species" differ- 
ing in local molecular structure: one of lower density and 
disorder, and the other of higher density and disorder 
Following the emergence of evidence for the continuum 
nature of the local structure and bonding in the liquid 
under ambient conditions, mixture models of water faded 
from prominence JJ. However, the proposal of a LLPT 
has renewed interest in mixture models, in which sponta- 
neous LDL- like and HDL-like structural fluctuations play 
the role of the mixed species [1, |^ . Recent experiments 
have been interpreted as evidence for such mixture-like 
fluctuations in water [To| , although this interpretation is 
disputed [TTj . 

One commonly discussed mixture model for water 
is based on the Gibbs free energy G of a binary regular 
solution, given by. 



G = Ga{1- X) + GbX + wX{l- X) 
+RT[X\tlX + [1 - X)\n{l -X)], 



(1) 



where T is the temperature, X is the concentration of 
component B, Ga and Gb are the respective free en- 
ergies of the pure A and pure B liquids, and R is the 
gas constant. The energy of mixing is quantified by 
w. It is further assumed that the two species can in- 
terconvert [A ^ B), and that the free energy difference 



between the two pure phases is given by Gb — Ga — 
AE - TAS + PAV. Here AE, AS and AV are respec- 
tively the difference in energy, entropy and molar volume 
of the two pure phases (e.g. AV = Vb — Va), and for 
simplicity are assumed to be constant with respect to T 
and pressure P. In this modified regular solution (MRS) 
model, the equilibrium value of X — x at fixed P and 
T is determined by minimizing Eq. [T] with respect to 
X. For w > 0, a liquid-liquid critical point occurs at 
Tc = w/2R, Pc = {TcAS - AE)/AV, and Xc = 1/2. 
This MRS model was originally developed by Rapaport 
to account for melting line maxima in pure liquids |12l |. 
and is a member of a large group of two-state models 
that have been applied to thermodynamic and relaxation 
phenomena in one-component liquids [l^ | and crys- 
tals [15„il6|]. 

MRS models have illustrated how a LLPT and the 
thermodynamic anomalies of water may be interrelated. 
These models have also been reported to be in quanti- 
tative agreement with experimental data for water d, |^, 

, but such comparisons are complicated by two signifi- 
cant challenges. First, the four model parameters w, AE, 
AS, and AV are not known for real water, and so they 
must be estimated indirectly, e.g. from the properties of 
the amorphous ices. Second, in order to compute thermo- 
dynamic properties from Eq. [1] the free energy function 
of a reference state [e.g. Ga{P, T)] is required. By them- 
selves, the four model parameters determine only the 
"anomalous" contribution to thermodynamic properties 
arising from the variation of x. Hence, to fit the model 
to experimental data without knowledge of Ga{P,T), a 
"normal" contribution must first be estimated and sub- 
tracted from the data. This process introduces additional 
and difficult-to-estimate parameters. Consequently, the 
regime of validity of MRS models for describing behav- 
ior near a LLPT, if such a regime exists at all, remains 
uncertain. 

In this Letter, we use simulation data to test if a ther- 
modynamic regime exists in which a MRS model can de- 
scribe a water-like liquid near a LLPT. We study the ST2 
model of water [l^ , as it provides a context in which the 
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FIG. 1: (a) T-P and (b) p-T projections of the properties of 
ST2 water. From Ref. we show the hne of density max- 
ima (thick black hne); density minima (thin black line); Kt 
maxima (thick blue line); Kt minima (thin blue line); the liq- 
uid spinodal (diamonds); the HDL spinodal (down triangles); 
and the LDL spinodal (up triangles). We also show the locus 
along which x — 0.5 (open red squares); the location of the 
critical point (open black circle); and the state points falling 
inside the region TZ defined in the text (crosses). The green 
line in (a) is a line of slope m, and in (b) is the coexistence 
curve predicted by the MRS model. 



two difBculties described above can be avoided. First, 
ST2 exhibits a well-characterized LLPT We can 

therefore determine the four model parameters of the 
MRS model directly. Second, we identify a property of 
the local structure in ST2 that allows us to estimate the 
concentrations of LDL-like and HDL-like species. We 
compare these concentrations directly to the predictions 
of the MRS model, thus avoiding the need to decompose 
thermodynamic properties into normal and anomalous 
contributions. 

Our results are based on molecular dynamics simula- 
tions of a system of = 1728 ST2 molecules. Our runs 
are conducted at fixed volume V, and T is controlled us- 
ing Berendsen's method [l9| . Long-range contributions 
to electrostatic interactions are approximated using the 
reaction field method. We study a wide range of states: 
from T ~ 220 to 400 K, in 5 K steps; and from density 
p — 0.8 to 1.1 g/cm'^, in steps of 0.01 g/cm^. Complete 
details of our simulation procedure are as described in 
Ref. J^. Fig.H] summarizes the known phase behavior 
of ST2. As reported in Ref. [l^l, a liquid-liquid critical 




FIG. 2: (a) gsir) at T = 245 K, for p = 0.88 to 1.06 g/cm^ 
in steps of 0.02 g/cm^. (b) 35 (r) at p = 0.96 g/cm^, for 
T = 220 to 270 K in steps of 10 K. For clarity, curves for 
T > 220 K are successively shifted upward by 0.05. (c) A 
snapshot of a system of TV = 13 824 ST2 molecules at T = 
245 K and p = 0.96 g/cm^; at this state we find P = 191 MPa. 
Blue molecules have rs > 0.35 nm; red molecules have rs < 
0.35 nm. 



point occurs in the vicinity of Tc = 245 K, Pc ~ 185 MPa 
and Pc = 0.94 g/cm'^. 

First, we seek a criterion for assigning molecules to 
LDL-like and HDL-like species. In Fig. [2j we analyze the 
liquid structure near the critical point in terms of the 
distance from the O atom of each molecule to its fifth- 
nearest neighbor. Following Ref. [2^, we define 55 (r) 
such that pg5{r) is the average density of fifth- nearest 
neighbors of an O atom at the origin, as found in a vol- 
ume element at a distance r. So defined, the conven- 
tional pair correlation function g(r) = X]i=i5n('")- 
the range of p studied here, fifth-nearest neighbors are 
located over a range of distances that span the first min- 
imum in the 0-0 pair correlation function, and thus 
is an indicator of the degree to which the tetrahedral 
structure of the first coordination shell is disrupted by 
additional neighbors. At T = Tc we find that rc, is typ- 
ically greater than 0.35 nm in the LDL phase (p < Pc), 
while in the HDL phase (p > Pc) is typically less than 
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FIG. 3: (a) Isotherms of function of P. (b) Isobars 

of X cLS cL function of T. Solid curves are the corresponding 
predictions of Eq. [T] 



0.35 nm [Fig. ^a)]. Near p = pc for T < T^, 55 (r) 
has a bimodal shape indicating the presence of distinct 
populations of LDL-Uke and HDL-hke coordination en- 
vironments [Fig. [2Ib)]. We therefore adopt as a local 
order parameter for assigning molecules to two species: 
"A" molecules are LDL-like and have rs > 0.35 nm; "B" 
molecules are HDL-like and have < 0.35 nm. Fig.[2fc) 
shows an equilibrium configuration from a separate sim- 
ulation of iV = 13 824 ST2 molecules at a state within 
error of the critical temperature and density. This im- 
age, in which A and B molecules are shown in different 
colors, confirms the presence of large, spatially correlated 
clusters of each species, as expected near a critical point. 

We evaluate x, the equilibrium number concentration 
of B molecules in the system, at each (p, T) state point as 
a time average over the equilibrium configurations gener- 
ated during the run. We also evaluate P, to allow us to 
analyze x{P,T) as well as x{p,T). Our results for x are 
shown as isotherms as a function of P in Fig. EJa), and 
as isobars as a function of T in Fig. ^h) . Approaching 
the critical point, the qualitative behavior of x is as ex- 
pected for a liquid mixture approaching a LLPT: Both 
isotherms (as P — Pc) and isobars (as T — > Tc) become 
more steeply sloped. 

Next we evaluate the parameters of the MRS model ap- 
propriate for ST2. We confirm the estimates of Rcf. [18] 
for Tc and Pc by examining two P-V isotherms straddling 



Tc (inset, Fig. |4]). Along the T = 250 K isotherm, P is 
monotonic in V, while the first sign of a "van der Waals 
loop" occurs along the T = 245 K isotherm. Given the 
scatter of the data along these isotherms, we estimate 
Tc = 247 ± 3 K and Pc = 185 ± 15 MPa. 

We estimate AV by noting that in the MRS model 
V = Vn + A.Vx, where Vn is the "normal" or non-singular 
contribution to V, and AVx is the "anomalous" contri- 
bution due to the variation of x [9]. Vn in general de- 
pends on both P and T. However, along the critical 
isotherm near Vc — ^/pc, both P and T are constant (in- 
set, Fig.S]), and therefore Vn is constant. Accordingly, we 
estimate AV from the slope of V versus x along the crit- 
ical isotherm in the interval 0.4 < a; < 0.6. As we will see 
below, this range of x spans the states near pc- By averag- 
ing the slopes (and their errors) obtained for T = 245 and 
250 K, we estimate AV = -5.0 ± 0.2 cmVmol (Fig. g]). 

As r ^ Tc from above, the MRS model predicts that 
the value of the isothermal compressibility Kt is increas- 
ingly dominated by a term containing {dx/dP)T, which 
diverges at the critical point [9]. In the model, the max- 
ima of isotherms of {dx/dP)T occur at x = 1/2 for 
all T. As a consequence, the locus along which Kt is 
a maximum, which converges to the "Widom line" as 
T ^ Tc [2l|, should also converge with the a; = 1/2 locus 
in the region approaching the critical point. The set of 
points satisfying a; = 1/2 is plotted in Figs.Ilfa) and (b) . 
along with the locus of Kt maxima reported in Ref. [l8| . 
We find that the x = 1/2 locus is in excellent agreement 
with the locus of Kt maxima for T < 290 K, and also 
that the location of the critical point is consistent with 
the model prediction of Xc = 1/2 in both the P-T and 
p-T planes. Computing the value of p at which a; = 1/2 
at T = Tc predicts pc = 0.955 ± 0.01 g/cm^. 

We also find that the a; = 1/2 locus, both above and 
below Tc, closely approximates a straight line in the P- 
T plane for T < 290 K. In the MRS model, both the 
X = 1/2 locus for T > Tc, as well as the coexistence 
curve for T < Tc, follow a straight line in the P-T plane 
given by P* = {AS/AV)T~{AE/AV) 0. From a linear 
fit to the a: = 1/2 locus from 250 to 280 K, we obtain the 
slope TO = AS/AV = -4.3 ± 0.2 MPa/K. We note that 
the resulting prediction for the coexistence line (with a 
Clapeyron slope of to) as expected lies between the LDL 
and HDL spinodal lines in Fig.[TJa). 

Our estimates for AV, Tc, Pc, and to completely de- 
termine the four model parameters {AV, w, AE, and 
AS) required for Eq. [1] We use these values to obtain 
the MRS model prediction for x{P,T), and compare the 
results to the data for x from our simulations (Fig. [3]). 
We find that inside a region TZ, defined approximately as 
0.25 <x< 0.75 and 250K < T < 290 K, the MRS model 
is in good agreement with the values of x computed from 
our simulations. At the boundaries of TZ and beyond, the 
agreement rapidly degrades. 

To test the robustness of the agreement between the 
MRS model and our data in the region TZ, we carry out 
a least-squares fit of the model to our data for x in this 
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FIG. 4: Isotherms near Tc of V as a function of x. The 
lines are linear fits to the data for 0.4 > a; > 0.6. For clar- 
ity, the data for T = 245 K have been shifted downward by 
0.5 cm^/mol. Inset: Isotherms of P versus V straddling Tc. 

region, allowing all four model parameters to vary. We 
select all distinct pairs of isotherms of a; in 7?. that are 
at least 15 K apart, and fit the model to each of the 21 
data subsets so defined. This gives 21 separate estimates 
for each fit parameter, from which we compute the mean 
and the standard deviation. The results are in excellent 
agreement with the values obtained above: Tc — 247 ± 
3 K, Pc = 181 ±11 MPa, AV^ = -5.2 ±0.6 cm^/mol, and 
m = -4.3 ±0.2 MPa/K. 

Our results thus demonstrate that a mixture model 
can indeed provide a quantitatively accurate description 
of a water-like liquid, in the specific case that the liquid 
exhibits a LLPT. Our MRS model successfully predicts 
the concentrations of LDL-like and HDL-like structural 
fluctuations in the region 72., which lies inside the lo- 
cus of density extrema, above Tc, and spans the range 
X = 0.5 ± 0.25 centered on the Widom line (Fig. [T|). The 
signatures for the onset of this "mixture-model regime" 
are the merging of the Widom line with the x = 1/2 lo- 
cus, and the observation of a linear Widom line in the 
P-T plane. These signatures may be useful for assessing 
other water-like liquids (either in simulations or experi- 
ments) for mixture-like behavior. 

Outside TZ, the MRS model does a poor job estimat- 
ing X. Fig. [IJa) shows that the high-T boundary of TZ 
is in the vicinity of 290 K, which is 90% of the high- 
est T (323 K) reached by the line of density maxima for 
ST2 water. For real water, 90% of the temperature of 
maximum density at ambient pressure (277 K) gives an 
estimate of 249 K (-24 C) for the highest T at which 
mixture-like behavior might be observed experimentally. 



This estimate supports the view that mixture models are 
not appropriate for interpreting the behavior of real wa- 
ter at ambient T ~ 300 K [ll|. 

We also note that the MRS model fails for T < T^. 
Fig. HKb) shows that the coexistence curve predicted by 
the MRS model [9] is too narrow since, in violation of 
thermodynamics, it lies inside the estimate of the LDL 
and HDL spinodal lines. The MRS model is a mean-field 
theory, and thus the shape of the coexistence curve near 
Tc obeys {p - p,) oc [{T, - T)lT,f, with /3 = 1/2. How- 
ever, the LLPT in ST2 water has been shown to belong 
to the 3D Ising universality class, for which /3 ~ 0.327 
Hence it is to be expected that the MRS model will un- 
derestimate the density difference between the coexisting 
phases as T decreases below Tc. In Fig. [3] we also note 
that there are significant deviations between the model 
and the data in the limits of large and small x for T > Tc, 
highlighting that the nearly pure A and B phases are 
poorly described by the model. This may also contribute 
to the discrepancy between the model and the data for 
the coexistence curve for T < Tc. 

However, within the mixture-model regime defined by 
72, our work shows that the examination of a local struc- 
tural property as a function of T and P can yield accu- 
rate information concerning the location of the Widom 
line and the critical point of a LLPT. In our case, the 
local structure is quantified in terms of a;, which is de- 
termined by the values of individual molecules. A 
number of simulation studies have used the behavior of 
55 (r) and related measures as evidence for a LLPT, in 
both tetrahedral 0, [2^ and non-tetrahedral liquids [l^l . 
Our results validate this approach, and further, confirm 
that the structures relevant to the LLPT in water-like 
liquids are highly localized, extending no farther than 
the second coordination shell. 

If a mixture-model regime exists for real water, our re- 
sults suggest that it will be found in the vicinity of the 
Widom line. States on the Widom line have yet to be 
studied in experiments on bulk supercooled water, due 
to the onset of rapid ice crystallization. However, for 
tetrahedral liquids in which the region of the Widom line 
is accessible, which may be the case for nanoconfined 
water [2^, our results demonstrate that probes of local 
molecular structure, in concert with mixture-model con- 
cepts, can be used to elucidate the properties of a LLPT 
and its associated critical point. 
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